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Abstract 

We  generalize  the  recently  introduced  twist- Kirchhoff  theory  of  rectangular  plate  ele¬ 
ments  to  arbitrary  quadrilateral  elements.  A  key  feature  is  the  use  of  Raviart-Thomas 
vector-field  approximations  for  rotations.  To  preserve  continuity  of  the  normal  compo¬ 
nents  of  the  rotation  vector  across  mesh  edges,  we  employ  the  Piola  transformation  to 
map  the  rotations  from  the  parent  domain  to  the  physical  domain.  These  elements  pos¬ 
sess  a  unique  combination  of  efficiency  and  robustness  in  that  minimal  quadrature  rules 
are  sufficient  to  guarantee  stability  without  rank  deficiency.  In  particular,  only  one-point 
Gauss  quadrature  is  required  for  the  lowest-order  element  in  the  twist-Kirchhoff  family. 
We  numerically  study  the  convergence  and  accuracy  of  the  first  two  members  of  the  twist- 
Kirchhoff  family  of  quadrilateral  elements  on  square,  rhombic  and  circular  plate  problems. 
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1  Introduction 


Considerable  research  effort  has  been  directed  toward  the  development  of  efficient  and 
reliable  quadrilateral  plate  finite  elements.  Quadrilateral  elements  are  particularly  at¬ 
tractive  for  the  discretization  of  plates  of  arbitrary  geometric  shapes. 

Most  studies  have  been  focused  on  the  Reissner-Mindlin  theory,  in  which  transverse 
shear  strains  are  accounted  for.  This  circumvents  difficulties  associated  with  C 1  require¬ 
ments  of  the  classical  Poisson-Kirchhoff  theory.  However,  standard  low-order  Reissner- 
Mindlin-type  elements  suffer  from  shear  locking  when  applied  to  thin  plates.  An  early 
remedy  to  this  was  the  reduced/selective  integration  technique  [12,  11],  Unfortunately, 
this  leads  to  rank  deficiency  and,  for  certain  boundary  conditions,  to  zero  energy  modes 
in  excess  of  the  three  physical  rigid  body  modes.  Many  research  efforts  have  been  under¬ 
taken  to  develop  Reissner-Mindlin  plate  elements  which  have  correct  rank  and  maintain 
accuracy  in  both  thin  and  thick  plate  applications,  see  e.g.  [13,  14,  4,  2,  3,  7,  6,  9]. 

Recently,  based  on  the  so-called  twist-Kirchhoff  theory,  which  conceptually  lies  “in 
between”  the  classical  Reissner-Mindlin  and  Poisson-Kirchhoff  theories,  a  new  family  of 
rectangular  plate  elements  has  been  proposed  for  the  analysis  of  thin  plates  [8].  A  salient 
feature  of  the  formulation  is  that  the  mixed  formulation  of  the  lowest-order  element 
for  the  limit  thin-plate  problem  is  exactly  integrated  with  one-point  Gauss  quadrature. 
The  element  has  no  rank  deficiency,  and  consequently  does  not  require  the  stabilization  of 
“hourglass”  or  any  other  singular  modes.  An  “equivalent”  primal  element,  eliminating  the 
transverse  shear  force  Lagrange  multipliers,  can  be  implemented  by  the  penalty  method. 
One-point  Gauss  quadrature  of  the  primal  element  attains  the  same  attributes  and,  in  the 
thin  plate  limit,  it  converges  to  the  mixed  element.  Consequently,  one-point,  quadrature 
is  also  viewed  as  exact  in  this  case.  This  element  attains  a  combination  of  stability 
and  efficiency  never  before  achieved  by  a  one-point  quadrature  quadrilateral  element. 
Consequently,  we  anticipate  that,  when  generalized  to  the  shell  setting,  it  may  offer 
significant  potential  for  economical  and  robust  explicit  crash  dynamics  and  sheet  metal 
forming  applications.  The  next  highest-order  element  shares  similar  properties.  Its  full 
integration  rule  is  2  x  2  Gauss  quadrature.  Likewise,  all  higher-order  elements  in  the 
family  behave  similarly. 

The  unique  mathematical  aspect  of  these  formulations  is  the  use  of  Raviart-Thomas 
vector-held  approximations  for  the  rotations.  We  employ  a  Piola.  transformation  to  de¬ 
fine  the  rotations  in  physical  space  for  the  arbitrary  quadrilateral  configuration.  This 
preserves  the  continuity  of  normal  rotations  across  element  interfaces.  We  generalize  the 
theory  to  account  for  the  Piola  transformation  and  show  how  to  invoke  the  twist-Kirchhoff 
hypothesis  in  the  arbitrary  quadrilateral  configuration.  Based  on  this,  we  implement  and 
numerically  study  the  convergence  behavior  of  the  first  two  elements  of  the  family  in 
distorted  configurations. 

It  is  well  known  that  the  approximation  properties  of  //(div)-conforming  Raviart- 
Thomas  elements  deteriorate  with  mesh  distortion  [1],  In  fact,  there  exist  sequences  of 
meshes  for  which  the  lowest-order  Raviart-Thomas  discretizations  of  the  mixed  Laplacian 
do  not  converge  with  respect  to  the  //(div)-norm.  This  problem  is  alleviated  through 
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the  utilization  of  higher-order  Raviart-Thomas  elements,  albeit  at  the  cost  of  a  power 
of  h.  This  being  said,  the  approximation  properties  of  lowest-order  Raviart-Thomas  el¬ 
ements  remain  optimal  for  regular  families  of  asymptotically  parallelogram  meshes  [5]. 
Furthermore,  it  is  manifestly  apparent  that  the  lowest-order  quadrilateral  element  may 
not  be  degenerated  to  a  triangle  in  the  twist- Kir chhoff  theory,  since  in  that  case  the 
transverse  displacement  field  becomes  linear  in  each  element  and  the  twist  component 
of  curvature,  which  is  calculated  as  the  mixed  second-derivative  of  the  transverse  dis¬ 
placement,  is  therefore  identically  zero.  Nevertheless,  we  investigate  the  behavior  of 
the  twist-Kirchhoff  elements  on  highly  distorted  configurations,  including  ones  in  which 
quadrilaterals  are  degenerated  to  triangles.  Indeed,  for  the  lowest-order  element,  mo¬ 
ments  in  degenerated  triangles  do  not  converge  to  correct  values.  However,  in  all  other 
cases,  the  lowest-order  element  behaves  well.  The  second-order  element  behaves  well  in 
all  cases,  including  the  case  of  degenerated  triangles. 

Square  plate  problems  are  discretized  with  asymptotically  parallelogram-shape-regular 
meshes  and  we  observe  the  same  optimal  rates  of  convergence  for  all  quantities  consid¬ 
ered  as  observed  for  rectangular  elements  in  our  previous  work.  Convergence  is  also 
noted  for  the  rhombic,  simply-supported  plate  problem,  known  as  Morley’s  problem, 
whose  solution  possesses  singularities  at  the  obtuse  vertices.  For  circular  plates  we  stud¬ 
ied  discretizations  in  which  singularities  were  introduced  in  the  geometrical  mappings  of 
elements.  In  every  case,  except  for  the  aforementioned  one,  convergence  was  obtained 
for  all  quantities  considered,  even  at  the  singularities  in  the  geometrical  mapping. 


2  Twist-Kirchhoff  Plate  Theory 

Let  H  C  1Z2  denote  the  mid-plane  of  a  plate  undergoing  rotation  6  =  [9X  9y]T  and 
transverse  deflection  w  as  a  result  of  applied  distributed  transverse  load  q  and  prescribed 
boundary  rotation  6  and  deflection  w.  Further,  let  dQ  =  r^rUr/)  denote  the  boundary  of 
f],  with  T n  and  T n  the  Neumann  and  Dirichlet  parts,  respectively,  such  that  Tn  fl  T ^  = 
0.  The  plate  is  assumed  to  be  linear  elastic,  homogeneous  and  isotropic  with  Young’s 
modulus  E ,  Poisson’s  ratio  v  and  constant  thickness  t. 

Let  us  also  introduce  the  kinematically  admissible  space  U  as 

U  =  {(0,  w)  E  -H1(h2)2  x  'H1(H)|  6  =  9  and  w  =  w  on  T#}  (1) 


In  the  classical  Reissner-Mindlin  formulation,  the  total  potential  energy  of  the  plate 


np :  U  — ^  72.  is  given  as  follows 


n p(0,w)=  I  ]-ktDmk,  +  ^7T£)y7  dtt  -  f 
Jn  *  *  Jn 


qw  dQ 


with  k  and  7  the  curvature  tensor  and  shear  strain  vector  of  the  plate  given  by 


^XX 

@x,x 

HL  = 

Kyy 

= 

®y,y 

1 

9x,y  +  9ytX 

and  7  = 


lx 

ly 


,x  @x 
w,y  ~ 


(2) 


(3) 
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and  Dm  and  Dy  the  bending/twisting  and  shear  constitutive  matrices  defined  as 


Dm  =  D 


1  1/  0 

v  10 
0  0 


and  Dy  =  kGt 


1  0 
0  1 


(4) 


where  D  =  Et3/(  12(1  —  z/2)),  G  =  E/( 2(1  +  u))  and  k  is  the  shear  correction  factor, 
usually  taken  to  be  5/6. 

The  minimization  of  the  total  potential  energy  (2)  leads  to  the  equilibrium  equations 
of  the  Reissner-Mindlin  plate  model. 

From  a  numerical  point  of  view,  the  difficulty  with  this  model  is  the  matching  of  the 
approximating  spaces  for  6  and  w.  As  t  tends  to  zero,  the  shear  strain  vector  7  must 
tend  to  zero  as  well,  i.e.,  the  shear  strains  =  w>x  —  9X  and  7^  =  wtV  —  9y  must  tend 
to  zero.  If  this  is  not  allowed  by  the  approximating  spaces,  the  result  is  a  deterioration 
of  the  numerical  solution,  well  known  in  the  literature  as  shear  locking.  The  situation  is 
particularly  vexing  if  we  wish  to  use  low-order  approximations. 

On  the  other  hand,  finite  elements  based  on  the  classical  Kirchhoff  theory  of  thin 
plate  bending  require  (^-continuity  of  the  transverse  displacement.  However,  whereas 
O0  finite  elements  of  various  shapes  and  numbers  of  nodes  are  readily  formulated,  the 
construction  of  multidimensional  C1  elements  has  proven  to  be  far  more  challenging. 

Recently,  a  new  theory,  the  twist- Kirchhoff  theory ,  has  been  proposed  for  thin  plates 
[8].  This  theory  lies  “in  between”  the  classical  Reissner-Mindlin  and  Poisson-Kirchhoff 
theories.  The  main  ingredient  of  this  theory  consists  in  replacing  the  Reissner-Mindlin 
twist  component  of  curvature,  defined  in  terms  of  first-order  derivatives  of  the  rotation 
components,  with  the  classical  Kirchhoff  twist  component  of  curvature,  defined  in  terms 
of  the  second-order  cross  derivative  of  the  transverse  displacement,  while  retaining  the 
definitions  of  the  bending  curvatures.  In  a  rectangular  Cartesian  frame,  this  amounts  to 
setting 

9 x,X 

K  =  9y}y 

,xy 

Based  on  this  theory,  a  new  family  of  rectangular  displacement/rotation-based  plate 
elements  has  been  proposed  for  the  analysis  of  thin  plates  [8].  This  new  family  of  elements 
uses  H (div)-conforming  Raviart- Thomas  vector  fields  of  order  r— 1  for  the  rotation  vector, 
where  r  >  1,  standard  C'°-continuous  piecewise  bi-Lagrange  functions  of  order  r  for  the 
transverse  displacement,  and  an  r  x  r-point  Gaussian  quadrature  rule  for  the  transverse 
shear  terms.  This  corresponds  exactly  to  a  mixed  formulation  in  which  the  transverse 
shear  force  resultants  are  assumed  to  be  discontinuous  piecewise  bi-Lagrange  functions  of 
order  r  —  1  over  each  element.  This  family  of  elements  has  been  mathematically  proven 
to  be  stable,  i.e.,  possessing  no  hourglass  modes,  and  to  be  free  from  shear-locking. 
Most  notably,  the  lowest-order  element  of  this  family  possesses  only  eight  degrees-of- 
freedom,  four  vertex  transverse  displacements  and  four  mid-side  rotations,  allowing  exact 
evaluation  of  the  stiffness  matrix  with  one-point  Gaussian  quadrature. 


(5) 
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We  present  in  this  paper  the  generalization  of  this  theory  to  the  arbitrary  quadrilateral 
element  case.  This  requires  special  treatment  of  the  Raviart-Thomas  rotation  field. 

Toward  this  end,  let  us  first  consider  a  nondegenerate  mapping  from  a  reference  square 
K  =  [—1, 1]  x  [—1, 1]  to  a  physical  domain  of  arbitrary  quadrilateral  shape  K  as  illustrated 
in  Figure  1,  where  </>  :  K  — >  K  is  a  regular  (orientation  preserving)  C 1  mapping,  and 
(x,  y)  and  (£,  rj)  denote  Cartesian  and  curvilinear  coordinates,  respectively. 


<t> 


Figure  1:  Geometric  mapping  from  the  parent  domain  to  the  physical  domain 


Further,  let  6  be  a  rotation  vector  field  on  K. 


The  Piola  transform  of  6  is  given  by 


G  =  -F0 

J 

where 

I 

On 

F  is  the  Jacobian  matrix  of  4>  defined  by 

F  =  [  x£  x,v 

_  y,t  y,v . 


r  r 

0  = 

_ i 

,  0  = 

(6) 

(7) 

(8) 


and  J  =  x^y>r)  —  xjT)y^  is  its  determinant.  The  Piola  transformation  preserves  continuity 
of  the  normal  components  of  the  rotation  vector  across  mesh  edges. 

The  transverse  deflection  is  mapped  using  a  standard  push-forward  with  the  relation¬ 
ship 

w  =  w  (9) 
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Making  use  of  these  relations,  the  classical  Kirchhoff  assumptions,  written  in  terms 
of  Cartesian  quantities  as 


u>,x  ~  Ox  =  0  (10a) 

W  y  ~  0y  =  0  (10b) 

can  be  rewritten  in  terms  of  curvilinear  quantities  as 

y,vw^  -  y,£W)V  -  x  -  x.ydr,  =0  (11a) 

-X,r,W,e  +  X£WlV  -  yjtz  -  y,rfi, 7  =0  (lib) 


where  we  have  cleared  the  factor  of  J”1. 

To  derive  the  equations  expressing  the  twist- Kirchhoff  hypothesis,  we  first  differentiate 
both  (11a)  and  (lib)  with  respect  to  £.  We  then  multiply  the  ^-derivative  of  (11a)  by  xjrt, 
and  the  ^-derivative  of  (lib)  by  ytV,  then  we  sum  and  solve  for  9V^.  The  result  is  (12a). 
We  proceed  in  analogous  fashion  to  obtain  (12b).  We  differentiate  (11a)  and  (lib)  with 
respect  to  rj,  then  we  multiply  the  former  result  by  and  the  latter  by  and  sum, 
and  solve  for  9^,  resulting  in  (12b). 


6^^  —  Aq[w  ^A\  +  JwM  —  9^A2  —  9^^A3  +  w  ^At±  —  9^A§^  (12a) 

%r,  =  B0(w:VB1  +  Jw^v  -  9t.B2  -  6v,vB3  +  w^B4  -  9VB5 )  (12b) 


where 


A)  =  1/  {x2v  +  y2v) 

-A  x^y^y  y^Xfy 

a2  x^x^y  ~f~  y-vy^v 

A  =  x£xiV  +  y,zy,v 
-A  y,ri%,££  %,ri  y,^ 

^5  x^x^  4-  y,vy,a 


and 


50  =  i/(®i  +  yi) 

B\  y^xxv  X£y,zn 

B2  '£,£•£, Zri  “i-  y,^y,iv 

B3  =  x^xiV  +  y.^y,v 
B*i  x^y  *) jjj  y  ^x  yy 

B§  x^x  yy  t  y,gym 


(13) 


Equations  (12a)  and  (12b)  represent  the  twist-Kirchhoff  assumptions  for  the  general 
quadrilateral  element  defined  in  the  curvilinear  system  (£,77).  See  Figure  1.  Note  that, 
in  the  rectangular  case,  (12a)  and  (12b)  simplify  to 


9V)i  =  (14a) 

%??  =  (14b) 


Further,  the  derivatives  of  the  Cartesian  components  of  the  rotation  vector  with 
respect  to  the  Cartesian  coordinates  can  be  expressed  in  terms  of  their  corresponding 
curvilinear  components  as  follows 


0x,x  =  J-2{{Cl-Dl)y,ri-{C2- 
9y,y  =  J-2((C4  -  DA)X£  ~  (C3  ~ 
9x,y  =  J-2{(C2-D2)x^-(C1- 
9y,x  =  J~2  {(C3  -  D3)ytV  -  (C4  - 


d2  )y,z) 

(15a) 

D3)xtV) 

(15b) 

Dl)X,y) 

(15c) 

d4  )y,() 

(15d) 
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with 


C\  —  +  x,t]9f]£  D\ 

( '2  (P (.m  ~L  '^,v^rhv  ^2 

C3  —  +  y,v@v,t.  -^3 

C4  y,tfii,rj  “F  U ■  tj^r/.r)  D4 


J  1(A±  —  Bi)Ei  —  Xfrdy  -  x^6j 
J  (R4  A  j )  E\  X^rjO^  X^pjOrq 
J~l(M  -  B1)E2  -  y.^nK  -  y.^Oj 
J  {^B-4  Ai)e2  y y,irtv^ri 


(16) 


where  ^  ^ 

E\  x  -b  x^Oyj  (17) 

E2  =  y,A  +  y,A 

Finally,  by  invoking  the  twist-Kirchhoff  assumptions  (12a)  and  (12b),  and  substituting 
them  into  (15a)-(15d),  we  obtain  the  components  of  the  twist-Kirchhoff  curvature  vector 
for  the  arbitrary  quadrilateral  element  case  as 


kxx  =  J  2 ((Ci  -  Djy#  -  (C2  -  D2)yA) 

Kyy  =  j~2((C4  ~  Da)x£  -  (C3  -  D3)xjV) 

KXy  —  J  2  ((C2  —  D2)x£  —  (Ci  —  Di)x>r)  +  (C3  —  D3)yv  —  (C4  —  D/j)y  ^j 
where  the  variables  Ci,  C2,  C3  and  C4  are  redefined  as 

Ci  =  x^9^^  +  x  ^Aq(w  ^A\  +  Jw^  —  9rjA2  —  9^^A3  +  w^A^  —  9^A^j 
C2  x  ^Bq  (w  ^Bi  T  Jw^  9^B2  9r/.ij  B  \  T  w  ^B^  T  x ^9^  ^ 

C3  =  ?/,£%£  +  y,rj^o(u>,^Ai  +  JwM  —  9r!A2  —  9^A3  +  w^A^  —  9^A^j 
C4  y  ^Bq{w  -qB\  T  Jw^  9^B2  9VtT1B3  T  w  ^B^  9^B 5)  A  y ^9^^ 


(18a) 

(18b) 

(18c) 


(19) 


3  Finite  Element  Approximations 

In  the  parent  domain,  the  rotations  are  represented  by  H (div)-conforming  Raviart- 
Tliomas  vector  fields  of  order  r  —  1,  where  r  >  1,  and  the  transverse  displacement  is 
represented  by  standard  C°-continuous,  piecewise  bi-Lagrange  functions  of  order  r.  Sim¬ 
ilarly  to  the  rectangular  element  case,  numerical  stability  in  the  arbitrary  quadrilateral 
element  case  requires  the  use  of  an  r  x  r-point  Gaussian  rule  on  the  transverse  shear 
term.  We  note  that,  in  the  rectangular  case,  the  r  x  r-point  Gaussian  rule  is  exact  for 
the  bending  term.  The  first  two  members  of  the  family  of  elements,  defined  in  the  parent 
domain,  are  schematically  illustrated  in  Figure  2. 

For  the  first-  and  second-order  elements,  he.,  for  r  =  1  and  r  =  2,  the  rotation 
components  and  the  transverse  displacement  fields  in  the  parametric  coordinates  are 
expressed  as 

2  6 


i=  1 

2 

i=l 

g 

i=l 

4 

and  = 

i— 1 

9 

(20) 

wh  =  J2  ^ Wi 

i=l 

wh  =  ^2  N™u>i 

i= 1 

7 


Figure  2:  First  two  members  of  the  element  family 

•  wh  degree  of  freedom 
=  0^  degree  of  freedom 

9^  degree  of  freedom 
X  Quadrature  points 
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respectively,  where  d\  and  9r‘  represent  mid-side  rotational  degrees-of- freedom  about  the 
r]-  and  ^-constant  coordinate  lines,  respectively,  whereas  wy  represent  vertex  displace¬ 
ment  degrees-of-freedom,  as  illustrated  in  Figure  2.  Note  that,  in  the  physical  domain, 
the  rotational  degrees-of-freedom  become  the  normal  components  of  the  rotation  on  the 
element  boundaries. 

In  the  first-order  element,  the  shape  functions  for  the  rotations  and  transverse  dis¬ 
placement  defined  in  parametric  coordinates  are 


<£  =  K1  -  f)  *r  =  i(i  -  0(i  -  o 

*?  =  1(1  +  0  JV?  =  j(i  +  0(1  —  0 

^"  =  1(1-0  ^3”  =  1(1 -0(1  +  0 

n£  =  1(1  +  0  %  =  i(i +  0(1  +  0 


whereas  in  the  second-order  element  they  are 


1){1  -  v) 
n9^  =  u  i-e)(i-V) 
N°3*  =  m+w-v) 

m+v) 
Nt*  =  l(i-e)(i+v) 
nS(  =  IZ(Z  +  W  +  v) 
N^  =  \v(  1-Ofo-l) 

fit"  =  k(  1-0(1-^2) 
n!v  =  H !- 0(^  +  1) 
fit"  =  &( !  + 0(^-1) 
St"  =  U  1  +  0(1-^2) 
Ntv  =  H  1  +  Ofo  +  l) 


AT  =  |^(1-  0(1  ~V) 

=  i-oa-o) 

%  =  -^(i  +  0(i-0 
i-m-v2) 
%  =  (  i-aa-^2) 

%  =  ie(i  +  0(i-T) 

%  =  -5^(l  “0(1 +  *7) 

%  =  ^(1+^(1 -a 

AT  =  i^(i  +  0(i  +  ??) 


(21) 


(22) 


4  Numerical  Results 


We  use  the  primal  formulation  with  a  reduced  quadrature  rule;  see  [8]  for  further  details. 
For  the  first-order  element  we  use  a  one-point  Gaussian  quadrature  rule,  whereas  for 
the  second-order  element  we  use  a  2  x  2-point  tensor-product  Gaussian  quadrature  rule. 
To  analyze  the  effectiveness  of  the  elements,  we  computed  solutions  for  the  first-  and 
second-order  cases  for  thicknesses  of  0.01,  0.001,  and  0.0001.  For  all  computations, 
a  shear  correction  factor  of  5/6  is  utilized  to  achieve  results  that  are  consistent  with 
classical  bending  theory  [10].  All  the  examples  refer  to  isotropic  plates  with  Young’s 
modulus  and  Poisson’s  ratio  taken  as  E  =  107  and  v  =  0.3.  A  uniform  distributed  load  of 
magnitude  q  =  1  is  applied  in  all  cases.  The  specification  of  essential  boundary  conditions 
is  identical  to  that  of  Poisson-Kirchhoff  theory,  he.,  for  simply-supported  edges  we  only 
impose  w  =  0,  for  clamped  edges  we  impose  w  =  0  and  9n  —  0.  The  geometrical  map  of 
the  element  utilizes  the  same  basis  functions  as  those  for  the  transverse  displacement. 
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The  energy  errors  presented  are  computed  with  the  bending  strain  energy  norm, 
defined  as 

||k||  =  J  k,tDmk  (23) 

As  proven  in  [8],  for  the  first-order  rectangular  element  we  have  second-order  conver¬ 
gence  for  w  in  L2(f2),  first-order  convergence  for  w  in  /A1  (0) ,  and  first-order  convergence 
for  |  At  I,  while  for  the  second-order  rectangular  element  we  have  third-order  convergence 
for  w  in  L2(Q),  second-order  convergence  for  w  in  Hl(Vt ),  and  second-order  convergence 
for  | /t | .  As  will  be  shown  next,  there  is  no  degradation  of  these  orders  of  convergence 
when  using  asymptotically  parallelogram-shape- regular  meshes,  be.,  meshes  whose  el¬ 
ements  converge  to  parallelograms  in  the  limit  of  mesh  refinement  and  whose  element 
aspect  ratios  are  uniformly  bounded  from  above  and  below. 

4.1  Square  Plate 

We  consider  an  isotropic  square  plate  subject  to  uniform  loading.  The  simply-supported 
case  is  presented  in  Section  4.1.1  and  the  fully-clamped  case  is  analyzed  in  Section  4.1.2. 

The  adopted  refinement  scheme  consists  of  “uniform”  refinements  of  a  non-uniform 
mesh  obtained  by  splitting  the  square  into  four  quadrilaterals,  as  illustrated  in  Figure  3, 
where  Dp  represents  the  distortion  parameter.  In  the  present  case,  Dp  was  set  to  1/10. 
Each  refinement  step  is  obtained  by  subdividing  each  quadrilateral  into  four  elements 
by  connecting  the  midpoints  of  the  opposite  edges.  This  leads  to  a  family  of  nei  x  nei 
asymptotically  parallelogram-shape- regular  meshes  as  shown  in  Figure  4  for  the  first 
three  meshes. 


Figure  3:  First  mesh  for  square  plate  problem  -  distortion  parameter 
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(a)  nei  =  2  (b)  nei  =  4  (c)  nei  =  8 

Figure  4:  Asymptotically  parallelogram  meshes  for  square  plate  problem 


4.1.1  Simply-Supported  Case 

Let  us  consider  the  simply-supported  square  plate  illustrated  in  Figure  5. 

V 


E  =  107 
v  =  0.3 
k  =  5/6 
L  =  1 

t  =  0.01,0.001,0.0001 
9  =  1 


x 

Figure  5:  Simply-supported  square  plate  model 


L 


L 


We  first  study  convergence  rates.  Since  there  is  no  known  analytical  solution  to  the 
twist-Kirchhoff  problem,  we  compare  discrete  solutions  with  a  heavily  refined  (128  x  128 
mesh)  solution  with  the  second-order  element.  We  use  this  to  determine  rates  of  con¬ 
vergence  in  integral  norms.  That  we  are  converging  to  the  exact  solution  will  be  verified 
by  tabular  results  presented  subsequently.  We  present  the  global  error  as  measured  by 
the  bending  strain  energy  norm  obtained  for  the  thickness  case  t  =  0.001  in  Figure  6. 
Similarly,  we  present  the  displacement  errors  as  measured  by  the  L2  and  H 1  norms  in 
Figures  7(a)  and  7(b),  respectively.  The  plotted  errors  are  normalized  by  the  norms  of 
the  reference  solution.  Examining  these  figures,  we  observe  rates  of  convergence  for  both 
the  first-  and  the  second-order  element  cases  which  are  the  same  as  the  optimal  rates  of 
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Figure  6:  The  normalized  total  error  produced  by  the  lowest-  and  second-order  plate 
elements  for  the  simply-supported  square  plate  under  a  uniform  load  for  a  thickness 
value  of  0.001. 


convergence  for  the  rectangular  element  case. 

We  now  study  the  convergence  of  the  center  displacement.  Tables  1  and  2  display 
the  convergence  of  the  center  displacement  for  the  first-  and  second-order  plate  elements 
respectively.  The  displayed  center  displacements  are  scaled  by  10 3D/  ( qL 4).  Comparing 
our  converged  twist-Kirchhoff  results  with  the  reference  Poisson-Kirchhoff  solution  [20], 
we  find  the  twist-Kirchhoff  center  displacement  converges  to  the  thin  plate  displacement 
from  above  as  the  thickness/width  ratio  t/L  — y  0,  as  we  might  have  anticipated.  To 
compare  our  twist-Kirchhoff  results  with  Reissner-Miudlin  theory,  we  have  simulated  a 
simply-supported  Reissner-Mindliu  plate  using  256  x  256  quadratic  Lagrange  rectangular 
elements  and  selective  reduced  integration  [11].  Table  3  displays  the  computed  center  dis¬ 
placements,  which  we  have  confirmed  are  converged  to  five  significant  digits.  We  find  the 
converged  twist-Kirchhoff  displacements  from  Tables  1  and  2  lie  below  the  corresponding 
Reissner-Mindlin  displacements  for  a  fixed  thickness/length  ratio  t/L.  This  result  also 
seems  consistent  with  the  “in  between”  nature  of  the  twist-Kichhoff  theory. 

We  next  study  the  convergence  of  the  center  bending  moment  about  the  x-axis.  Since 
the  discrete  center  bending  moment  is  not  well-defined,  we  sample  it  at  a  quadrature  point 
lying  closest  to  the  center  of  the  plate.  Tables  4  and  5  display  the  convergence  of  the  center 
moments  for  the  first-  and  second-order  plate  elements,  respectively.  The  displayed  center 
moments  are  scaled  by  102/  ( qL 2).  As  expected,  the  convergence  rate  for  the  bending 
moment  is  slower  than  that  for  the  center  displacement.  This  is  significantly  influenced 
by  the  fact  that  the  location  at  which  we  sample  the  bending  moment  within  an  element 
is  O(h),  where  h  &  l/nei,  and  thus  we  would  have  no  reason  to  expect  convergence  to 
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Figure  7:  (a)  The  normalized  L2  norm  of  the  displacement  produced  by  the  lowest-  and 
second-order  plate  elements  for  the  simply-supported  square  plate  under  a  uniform  load 
for  a  thickness  value  of  0.001.  (b)  The  normalized  H 1  norm  of  the  displacement  for  the 
same  problem. 


13 


Mesh 

t/L  =  0.01 

o 

o 

o 

I—1 

t/L  =  0.0001 

2x2 

3.86286 

3.86128 

3.86126 

4x4 

4.08572 

4.08479 

4.08479 

8x8 

4.07035 

4.06952 

4.06951 

16  x  16 

4.06500 

4.06416 

4.06415 

32  x  32 

4.06366 

4.06281 

4.06280 

64  x  64 

4.06332 

4.06247 

4.06247 

Table  1:  Center  displacement  (w  x  103  D  /  (qL4))  for  first-order  element  simply-supported 
square  plate  solutions  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  4.06235  [20]. 


Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

4.20035 

4.19957 

4.19956 

4x4 

4.07131 

4.07050 

4.07049 

8x8 

4.06369 

4.06286 

4.06285 

16  x  16 

4.06323 

4.06239 

4.06238 

32  x  32 

4.06321 

4.06236 

4.06235 

64  x  64 

4.06321 

4.06236 

4.06235 

Table  2:  Center  displacement  ( w  x  10 3D/(qL4))  for  second-order  element  simply- 
supported  square  plate  solutions  for  the  three  thickness/lcngth  ratios.  Reference  thin 
plate  limit  solution  is  4.06235  [20]. 


t/L  =  0.01  t/L  =  0.001  t/L  =  0.0001 
4.09930  4.06585  4.06268 


Table  3:  Center  displacement  (w  x  10 3D/(qL4))  for  simply-supported  Reissner-Mindlin 
plate.  Computed  using  256  x  256  quadratic  Lagrange  rectangular  elements  and  selective 
reduced  integration  [11], 
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Mesh 

o 

o 

o 

o 

o 

I—1 

t/L  =  0.0001 

2x2 

1.88468 

1.88516 

1.88517 

4x4 

4.17328 

4.17385 

4.17385 

8x8 

4.64061 

4.64103 

4.64104 

16  x  16 

4.75203 

4.75238 

4.75238 

32  x  32 

4.77916 

4.77969 

4.77968 

64  x  64 

4.78584 

4.78643 

4.78641 

Table  4:  “Center”  bending  moment  about  the  x-axis  (Mx  x  102/(gL2))  for  first-order 
element  simply-supported  square  plate  solutions  for  the  three  thickness/length  ratios. 
Reference  thin  plate  limit  solution  is  4.78864  [20]. 


Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

4.52602 

4.52648 

4.52649 

4x4 

4.68831 

4.68874 

4.68875 

8x8 

4.76160 

4.76167 

4.76167 

16  x  16 

4.78166 

4.78171 

4.78170 

32  x  32 

4.78646 

4.78694 

4.78688 

64  x  64 

4.78765 

4.78823 

4.78820 

Table  5:  “Center”  bending  moment  about  the  x-axis  (Mx  x  102/(gL2))  for  second-order 
element  simply-supported  square  plate  solutions  for  the  three  thickness/length  ratios. 
Reference  thin  plate  limit  solution  is  4.78864  [20]. 


be  any  better  than  first-order.  However,  we  observe  second-order  convergence  for  both 
elements,  with  the  absolute  values  of  the  error  being  more  accurate  for  the  second-order 
element. 

4.1.2  Fully-Clamped  Case 

Let  us  now  consider  the  case  of  a  fully-clamped  square  plate,  as  illustrated  in  Figure  8. 

In  Figure  9  we  present  the  global  error  as  measured  by  the  bending  strain  energy 
norm.  Similarly,  we  present  the  displacement  errors  as  measured  by  the  L2  and  H 1 
norms  in  Figures  10(a)  and  10(b),  respectively.  The  plotted  errors  are  normalized  by  the 
norms  of  the  exact  solution.  As  in  the  simply-supported  plate  problem,  we  observe  no 
degradation  of  the  optimal  rates  of  convergence  for  the  rectangular  element  case. 

Tables  6  and  7  display  the  convergence  of  the  center  displacement  for  the  first-  and 
second-order  plate  elements,  respectively.  As  can  be  seen,  the  twist-Kirchhoff  center  dis¬ 
placement  converges  to  the  thin  plate  displacement  from  above  as  the  thickness/length 
ratio  t/L  0.  To  compare  our  twist-Kirchhoff  results  with  Reissner-Mindlin  theory,  we 
simulated  a  clamped  Reissner-Mindlin  plate  using  256  x  256  quadratic  Lagrange  rect- 
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y 


E  =  107 
v  =  0.3 
k  =  5/6 
L  =  1 

t  =  0.01,0.001,0.0001 
9  =  1 


x 

Figure  8:  Fully-clamped  square  plate  model 


Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

0.00343 

0.00003 

0.00000 

4x4 

1.21557 

1.21318 

1.21316 

8x8 

1.26333 

1.26139 

1.26137 

16  x  16 

1.26653 

1.26471 

1.26469 

32  x  32 

1.26700 

1.26520 

1.26518 

64  x  64 

1.26710 

1.26530 

1.26529 

Table  6:  Center  displacement  (w  x  10 3D/(qL4))  for  first-order  element  fully-clamped 
square  plate  solutions  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  1.26532  [19]. 


angular  elements  and  selective  reduced  integration  [11],  Table  8  displays  the  computed 
center  displacements,  which  we  have  confirmed  are  converged  to  six  significant  digits. 
As  was  the  case  for  the  simply-supported  plate,  we  find  the  converged  twist-Kirchhoff 
displacements  lie  below  the  corresponding  Reissner-Mindlin  displacements  for  a  fixed 
thickness/length  ratio  t/L.  This  result  again  seems  consistent  with  the  “in  between” 
nature  of  the  twist- Kichhoff  theory. 

We  next  study  the  convergence  of  the  center  bending  moment  about  the  x-axis.  As 
was  done  for  the  simply-supported  case,  we  sample  the  discrete  bending  moment  at 
a  quadrature  point  lying  closest  to  the  center  of  the  plate.  Tables  9  and  10  display 
the  convergence  of  the  center  moment  for  the  first-  and  second-order  plate  elements, 
respectively.  The  displayed  center  moments  are  scaled  by  102/  ( qL 2).  We  note  that  the 
bending  moment  exhibits  second-order  convergence  for  both  the  first-  and  second-order 
case,  but  with  the  second-order  case  being  more  accurate  on  an  absolute  basis. 
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Figure  9:  The  normalized  total  error  produced  by  the  lowest-  and  second-order  plate 
elements  for  the  fully-clamped  isotropic  plate  under  a  uniform  load  for  a  thickness  value 
of  0.001. 


Mesh 

OR 

o 

o 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

1.50167 

1.49994 

1.49992 

4x4 

1.27850 

1.27664 

1.27662 

8x8 

1.26785 

1.26603 

1.26601 

16  x  16 

1.26718 

1.26538 

1.26536 

32  x  32 

1.26714 

1.26534 

1.26532 

64  x  64 

1.26714 

1.26534 

1.26532 

Table  7:  Center  displacement  (w  x  10 3D/(qL4))  for  second-order  element  fully-clamped 
square  plate  solutions  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  1.26532  [19]. 


t/L  =  0.01  t/L  =  0.001  t/L  =  0.0001 
1.26787  1.26534  1.26532 


Table  8:  Center  displacement  (w  x  10 3D/(qL4))  for  clamped  Reissner-Mindlin  plate. 
Computed  using  256  x  256  quadratic  Lagrange  rectangular  elements  and  selective  reduced 
integration  [11], 
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Figure  10:  (a)  The  normalized  L2  norm  of  the  displacement  produced  by  the  lowest-  and 
second-order  plate  elements  for  the  fully-clamped  isotropic  plate  under  a  uniform  load 
for  a  thickness  value  of  0.001.  (b)  The  normalized  H 1  norm  of  the  displacement  for  the 
same  problem. 
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Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

-0.00019 

0.00000 

0.00000 

4x4 

2.05023 

2.05186 

2.05188 

8x8 

2.22775 

2.23134 

2.23138 

16  x  16 

2.27135 

2.27615 

2.27630 

32  x  32 

2.28529 

2.28656 

2.28701 

64  x  64 

2.28903 

2.28928 

2.28962 

Table  9:  “Center”  bending  moment  about  the  x-axis  (Mx  x  102/ ( qL 2))  for  first-order  ele¬ 
ment  fully-clamped  square  plate  solutions  for  the  three  thickness/length  ratios.  Reference 
thin  plate  limit  solution  is  2.29051  [19]. 


Mesh 

oh 

o 

o 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

1.96176 

1.96023 

1.96022 

4x4 

2.16743 

2.16461 

2.16458 

8x8 

2.26193 

2.25608 

2.25595 

16  x  16 

2.28422 

2.28206 

2.28161 

32  x  32 

2.28879 

2.28881 

2.28828 

64  x  64 

2.28990 

2.29013 

2.29000 

Table  10:  “Center”  bending  moment  about  the  rc-axis  (Mx  x  102/(gL2))  for  second- 
order  element  fully-clamped  square  plate  solutions  for  the  three  thickness/lcngth  ratios. 
Reference  thin  plate  limit  solution  is  2.29051  [19]. 
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4.2  Simply-Supported  Morley’s  Plate 

Let  us  now  consider  an  isotropic  simply-supported  rhombic  plate  under  uniform  trans¬ 
verse  loading  as  illustrated  in  Figure  11.  In  this  case,  a  singular  behaviour  is  known  to 
occur  at  the  obtuse  vertices.  As  a  result,  we  anticipate  some  difficulties  may  arise.  In 
fact,  some  thin  plate  elements  yield  pathological  results  for  this  problem  [18].  An  ana¬ 
lytical  solution  given  as  an  infinite  series  was  obtained  in  [15]  for  the  Poisson-Kirchhoff 
case. 

E  =  107 
v  =  0.3 
k  =  5/6 
L  =  1 

t  =  0.01,0.001,0.0001 

q  =  i 

Figure  11:  Simply-supported  Morley’s  plate  model 

The  adopted  refinement  scheme  consists  in  uniform  nei  x  nei  parallelogram-shape- 
regular  meshes  as  shown  in  Figure  12  for  the  first  three  meshes. 

We  first  study  the  convergence  of  the  center  displacement.  Tables  11  and  12  display 
the  convergence  of  the  center  displacement  for  the  first-  and  second-order  plate  elements, 
respectively.  Comparing  our  converged  twist-Kirchhoff  results  with  the  reference  Poisson- 
Kirchhoff  solution  [15]  evaluated  retaining  nine  terms  of  the  infinite  series,  we  find  the 
twist-Kirchhoff  center  displacement  converges  to  the  thin  plate  displacement  from  above 
as  the  thickness/length  ratio  t/L  — y  0. 

We  next  study  the  convergence  of  the  maximum  and  minimum  principal  center  bend¬ 
ing  moments.  We  sample  the  discrete  bending  moment  at  a  quadrature  point  lying  closest 
to  the  center  of  the  plate.  Tables  13  and  14  display  the  convergence  of  the  center  bend¬ 
ing  moments  for  the  first-  and  second-order  elements,  respectively.  The  displayed  center 
moments  are  scaled  by  102/  ( qL 2).  We  note  that  the  bending  moment  converges  slowly 
for  both  the  first-  and  second-order  cases.  The  second-order  discretization  is  slightly 


(a)  nei  =  2  (b)  nei  =  4  (c)  nei  =  8 

Figure  12:  Parallelogram  meshes  for  Morley’s  plate  problem 
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Mesh 

t/L  =  0.01  t/L  =  0.001 

t/L  =  0.0001 

2x2 

0.26392 

0.26307 

0.26306 

4x4 

0.43539 

0.43473 

0.43472 

8x8 

0.43717 

0.43675 

0.43674 

16  x  16 

0.42091 

0.41935 

0.41933 

32  x  32 

0.41810 

0.41325 

0.41307 

64  x  64 

0.41794 

0.41177 

0.41066 

128  x  128 

0.41790 

0.41170 

0.40958 

Table  11:  Center  displacement  (w  x  10 3D/(qL4))  for  first-order  element  simply-suported 
Mor ley’s  plate  solutions  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  evaluated  retaining  nine  terms  in  the  series  is  0.40777  [15]. 


Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

0.50100 

0.50035 

0.50034 

4x4 

0.44586 

0.44545 

0.44544 

8x8 

0.42174 

0.42064 

0.42063 

16  x  16 

0.41799 

0.41379 

0.41367 

32  x  32 

0.41790 

0.41182 

0.41098 

64  x  64 

0.41789 

0.41169 

0.40972 

128  x  128 

0.41789 

0.41169 

0.40939 

Table  12:  Center  displacement  ( w  x  10 3H/(gL4))  for  second-order  element  simply- 
supported  Morley’s  plate  solutions  for  the  three  thickness/ length  ratios.  Reference  thin 
plate  limit  solution  evaluated  retaining  nine  terms  in  the  series  is  0.40777  [15]. 


more  accurate  in  absolute  terms.  These  low  rates  of  convergence  can  be  attributed  to 
the  singularities  at  the  obtuse  corners,  and  the  fact  that  the  meshes  arc  uniform  and  not 
graded  to  better  represent  the  singular  behavior. 

4.3  Circular  Plate 

Let  us  now  consider  an  isotropic  circular  plate  under  uniform  transverse  loading.  As  in 
the  square  plate  problem,  we  analyze  both  the  simply-supported  and  the  fully-clamped 
cases;  see  Sections  4.3.1  and  4.3.2,  respectively. 

Two  different  mesh  refinement  schemes  are  adopted  in  this  case.  One  is  based  on  a 
sequence  of  polar  meshes  (PM),  whereas  the  other  is  based  on  a  sequence  of  meshes  ob¬ 
tained  by  mapping  a  square  domain  onto  a  circle  (SCM).  We  note  that  for  the  PM-meshes, 
we  can  interpret  the  triangular  elements  around  the  center  as  degenerate  quadrilaterals 
with  one  edge  collapsed.  Hence,  from  the  computational  point  of  view,  for  the  coinci¬ 
dent  nodes,  the  associated  degrees-of- freedom  were  made  equivalent,  while  the  rotational 
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t/L  = 

0.01 

t/L  = 

0.001 

t/L  = 

0.0001 

Mesh 

Mmax 

M-min 

M-max 

M-min 

M-max 

Mmin 

2x2 

0.80932 

0.19377 

0.80935 

0.19410 

0.80935 

0.19410 

4x4 

1.67762 

0.47686 

1.67774 

0.47685 

1.67774 

0.47685 

8x8 

1.93080 

0.76694 

1.93162 

0.76060 

1.93163 

0.76054 

16  x  16 

1.94468 

0.97043 

1.95618 

0.90315 

1.95637 

0.90218 

32  x  32 

1.93018 

1.10883 

1.94570 

0.97020 

1.94832 

0.95879 

64  x  64 

1.93288 

1.11883 

1.92356 

1.06029 

1.93762 

0.99155 

128  x  128 

1.93374 

1.12036 

1.91706 

1.09950 

1.92559 

1.02745 

Table  13:  Maximum  and  minimum  principal  “center”  bending  moments  (iff  x  102/(gL2)) 
for  first-order  element  simply-supported  Morley’s  plate  solutions  for  the  three  thick¬ 
ness/length  ratios.  Reference  thin  plate  limit  solutions  evaluated  retaining  nine  terms  in 
the  series  are  1.9080  and  1.0834  [15]. 


Figure  13:  PM  meshes  with  bi-linear  mapping  for  circular  plate  problem 


degrees-of- freedom  associated  with  collapsed  edges  were  removed  from  the  system,  i.e., 
set  to  zero.  We  are  particularly  interested  in  the  effect  of  degenerating  the  elements  on 
local  derivative  quantities,  especially  the  bending  moments. 

For  the  PM-meshes  we  exploit  symmetry  conditions  and  only  discretize  the  first  quad¬ 
rant  of  the  plate,  whereas  for  the  SCM-meshes  we  discretize  the  whole  domain,  see  Figures 
13  and  14,  respectively.  nei  represents  in  this  case  the  adopted  number  of  elements  per 
side  in  one  quadrant  of  the  plate.  For  the  meshes  with  second-order  elements  we  use 
a  geometrical  bi-quadratic  mapping,  in  which  initially  straight  edges  in  the  reference 
domain  are  mapped  onto  curved  edges  in  the  physical  domain;  see  Figures  15  and  16. 
As  mentioned  in  Section  1,  we  have  no  illusions  about  the  lowest-order  twist-Kirchhoff 
element  being  able  to  deal  with  polar  meshes. 
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(a)  nei  =  2  (b)  nei  =  4  (c)  nei  =  8 

Figure  14:  SCM  meshes  with  bi-linear  mapping  for  circular  plate  problem 


Figure  15:  PM  meshes  with  bi-quadratic  mapping  for  circular  plate  problem 


(a)  nei  =  2  (b)  nez  =  4  (c)  ne;  =  8 

Figure  16:  SCM  meshes  with  bi-quadratic  mapping  for  circular  plate  problem 
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t/L  = 

0.01 

t/L  = 

0.001 

t/L  = 

0.0001 

Mesh 

Mmax 

M-min 

M-max 

M-min 

M-max 

■Mmin 

2x2 

2.00063 

1.27113 

2.00090 

1.27079 

2.00090 

1.27079 

4x4 

1.88260 

1.40727 

1.88092 

1.41116 

1.88090 

1.41119 

8x8 

1.88789 

1.28405 

1.87024 

1.32483 

1.87001 

1.32534 

16  x  16 

1.92966 

1.13306 

1.88602 

1.24460 

1.88357 

1.25070 

32  x  32 

1.93351 

1.12083 

1.90501 

1.15871 

1.88817 

1.20602 

64  x  64 

1.93390 

1.12085 

1.91713 

1.10084 

1.89412 

1.16676 

128  x  128 

1.93400 

1.12086 

1.91731 

1.10010 

1.90766 

1.10805 

Table  14:  Maximum  and  minimum  principal  “center”  bending  moments  (iff  x  102/(gL2)) 
for  second-order  element  simply-supported  Morley’s  plate  solutions  for  the  three  thick¬ 
ness/length  ratios.  Reference  thin  plate  limit  solutions  evaluated  retaining  nine  terms  in 
the  series  are  1.9080  and  1.0834  [15]. 

4.3.1  Simply-Supported  Case 

We  consider  a  simply-supported  circular  plate  of  radius  R  =  1,  as  illustrated  in  Figure 
17. 

For  the  PM-meshes,  the  only  essential  boundary  condition  imposed  on  the  curved 
portion  of  the  boundary  is  w  =  0,  while  on  the  vertical  and  horizontal  edges  of  the 
quadrant  boundary  the  only  essential  boundary  condition  imposed  is  9n  =  0.  For  the 
SCM-meshes,  the  essential  boundary  condition  imposed  on  the  whole  boundary  is  w  =  0. 


E  =  107 
v  =  0.3 
k  =  5/6 
R  =  1 

t  =  0.01,0.001,0.0001 
q  =  1 


Figure  17:  Simply-supported  circular  plate  model 

We  are  particularly  interested  in  studying  the  convergence  of  the  central  deflection, 
central  bending  moment  about  the  x-axis,  and  bending  moment  about  the  y-axis  at 
point  (x,  y )  =  (R,  0)  of  the  boundary.  The  exact  Poisson-Kirchhoff  solutions  for  these 
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t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.82936 

0.94349 

0.98339 

0.99542 

0.99882 

0.99976 


0.81177 

0.95153 

0.98787 

0.99704 

0.99933 

0.99991 


0.82923 

0.94337 

0.98327 

0.99531 

0.99871 

0.99965 


0.81166 

0.95143 

0.98777 

0.99694 

0.99924 

0.99981 


0.82923 

0.94337 

0.98327 

0.99531 

0.99871 

0.99966 


0.81166 

0.95143 

0.98777 

0.99694 

0.99923 

0.99981 


Table  15:  Center  displacement  (w  x  64(1  +  v)D/{{ 5  +  v)qR 4))  for  first-order  element 
simply-suported  circular  plate  solutions  for  the  three  thickness/length  ratios.  Reference 
thin  plate  limit  solution  is  1.0  [17]. 


t/R  = 

0.01 

t/R  = 

0.001 

t/R  = 

0.0001 

Mesh 

PM 

SCM 

PM 

SCM 

PM 

SCM 

2x2 

0.99441 

0.99947 

0.99423 

0.99937 

0.99422 

0.99937 

4x4 

0.99982 

1.00006 

0.99967 

0.99996 

0.99967 

0.99995 

8x8 

1.00010 

1.00010 

0.99998 

1.00000 

0.99998 

1.00000 

16  x  16 

1.00011 

1.00010 

1.00000 

1.00000 

0.99999 

1.00000 

32  x  32 

1.00011 

1.00010 

1.00000 

1.00000 

0.99992 

1.00000 

64  x  64 

1.00011 

1.00010 

1.00000 

1.00000 

0.99992 

1.00000 

Table  16:  Center  displacement  (w  x  64(1  +  u)D/{{ 5  +  u)qR1))  for  second-order  element 
simply-supported  circular  plate  solutions  with  bi-quadratic  mapping  for  the  three  thick¬ 
ness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 

quantities  are  given  by  [17]: 


fb  +  v\ 

ro(0-°)=(1  +  1,J64  D 

(24a) 

Mx( 0,  0)  —  (3  +  v) 

16 

(24b) 

MV(R,  0)  =  (1  — 

(24c) 

The  convergence  of  the  normalized  center  deflection  obtained  with  the  first-  and 
second-order  elements  is  displayed  in  Tables  15  and  16,  respectively.  Normalization  is 
taken  with  respect  to  the  exact  Poisson-Kirchhoff  solution  (24a). 

We  next  study  the  convergence  of  the  normalized  center  bending  moment  about  the 
x-axis.  We  sample  the  discrete  bending  moment  at  a  quadrature  point  lying  closest  to 
the  center  of  the  plate.  Tables  17  and  18  display  the  convergence  of  the  center  bending 
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t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.41008 

0.34687 

0.32177 

0.31372 

0.31129 

0.31058 


0.75691 

0.93594 

0.98378 

0.99593 

0.99898 

0.99974 


0.41008 

0.34687 

0.32177 

0.31372 

0.31129 

0.31058 


0.75690 

0.93594 

0.98378 

0.99593 

0.99898 

0.99975 


0.41008 

0.34687 

0.32177 

0.31372 

0.31129 

0.31058 


0.75690 

0.93594 

0.98378 

0.99593 

0.99898 

0.99974 


Table  17:  “Center”  bending  moment  about  the  x-axis  (Mx  x  16/(gi?2(3  +  is)))  for  first- 
order  element  simply-supported  plate  solutions  for  the  three  thickness/length  ratios.  Ref¬ 
erence  thin  plate  limit  solution  is  1.0  [17]. 


moment  obtained  using  the  two  refinement  schemes  for  the  first-  and  second-order  plate 
elements,  respectively.  The  displayed  center  moments  are  scaled  by  the  exact  Poisson- 
Kirchhoff  solution,  given  in  this  case  by  (24b).  We  note  that  the  bending  moment 
converges  rapidly  for  the  second-order  discretization.  However,  for  the  polar  mesh  case, 
the  first-order  element  bending  moments  do  not  converge  to  the  true  solutions.  This  may 
be  attributed  to  the  anticipated  difficulties  previously  alluded  to.  To  verify  that  this  is 
a  problem  attributed  to  the  lowest-order  Raviart- Thomas  vector  fields,  we  solved  the 
simply-supported  circular  plate  problem  with  the  standard  selectively  integrated  bilinear 
Reissner-Mindlin  element  and  found  that  the  moments  did  in  fact  converge.  This  issue 
also  occurs  for  the  fully-clamped  case;  see  the  next  section. 

To  conclude  the  analysis  of  the  present  example,  we  study  the  convergence  of  the 
bending  moment  about  the  y-axis  at  point  (x,  y)  =  (R,  0)  of  the  boundary.  Tables  19 
and  20  display  the  convergence  of  the  bending  moment  obtained  using  the  two  refine¬ 
ment  schemes  for  the  first-  and  second-order  plate  elements,  respectively.  The  displayed 
moments  are  scaled  by  (24c).  Examining  Table  19,  we  observe  that,  unlike  the  central 
bending  moments  obtained  with  the  polar  discretization  scheme,  no  convergence  issues 
seem  to  occur  for  the  bending  moment  about  the  y-axis  at  point  (x,  y)  =  (R,  0)  of  the 
boundary. 

4.3.2  Fully-Clamped  Case 

Let  us  now  consider  the  fully-clamped  circular  plate  illustrated  in  Figure  18. 

For  the  PM-meshes,  essential  boundary  conditions  w  =  9n  =  0  are  imposed  on  the 
curved  portion  of  the  boundary,  while  on  the  vertical  and  horizontal  edges  of  the  quadrant 
boundary  the  only  essential  boundary  conditions  imposed  are  9n  =  0.  For  the  SCM- 
meshes,  essential  boundary  conditions  w  =  9n  =  0  are  imposed  on  the  whole  boundary. 

We  are  interested  in  this  case  in  the  convergence  of  the  central  deflection,  central 
bending  moment  about  the  x-axis,  and  bending  moment  about  the  x-axis  at  point  (x,  y)  = 
(R,  0)  of  the  clamped  edge.  The  exact  Poisson-Kirchhoff  solutions  for  these  quantities 
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t/R  = 

0.01 

t/R  = 

0.001 

t/R  = 

0.0001 

Mesh 

PM 

SCM 

PM 

SCM 

PM 

SCM 

2x2 

0.88524 

0.98033 

0.88251 

0.98033 

0.88247 

0.98033 

4x4 

0.97411 

0.99547 

0.97026 

0.99547 

0.97011 

0.99547 

8x8 

0.99379 

0.99889 

0.99296 

0.99889 

0.99250 

0.99889 

16  x  16 

0.99845 

0.99972 

0.99843 

0.99972 

0.99814 

0.99972 

32  x  32 

0.99961 

0.99993 

0.99961 

0.99993 

0.99938 

0.99993 

64  x  64 

0.99990 

0.99998 

0.99990 

0.99998 

0.99970 

0.99998 

Table  18:  “Center”  bending  moment  about  the  x-axis  (Mx  x  16/ (qR2(3  +  is)))  for  second- 
order  element  simply-supported  plate  solutions  with  bi-quadratic  mapping  for  the  three 
thickness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


1.20381 

1.19732 

1.13215 

1.07526 

1.04000 

1.02060 


1.32128 

1.27713 

1.15872 

1.08262 

1.04179 

1.02090 


1.20381 

1.19732 

1.13215 

1.07526 

1.04000 

1.02060 


1.32136 

1.27727 

1.15890 

1.08282 

1.04200 

1.02111 


1.20381 

1.19732 

1.13215 

1.07526 

1.04000 

1.02061 


1.32136 

1.27727 

1.15890 

1.08282 

1.04200 

1.02112 


Table  19:  Bending  moment  about  the  y- axis  at  point  (x,  y )  =  (R,  0)  of  the  boundary 
(My  x  8/(qR2(l  —  u)))  for  first-order  element  simply-supported  plate  solutions  for  the 
three  thickness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


1.22662 

1.13080 

1.06880 

1.03516 

1.01775 

1.00892 


1.28616 

1.14230 

1.07086 

1.03540 

1.01763 

1.00872 


1.22619 

1.13062 

1.06879 

1.03515 

1.01775 

1.00892 


1.28631 

1.14247 

1.07106 

1.03561 

1.01785 

1.00894 


1.22618 

1.13061 

1.06878 

1.03515 

1.01769 

1.00886 


1.28631 

1.14247 

1.07106 

1.03562 

1.01786 

1.00894 


Table  20:  Bending  moment  about  the  y- axis  at  point  (x,  y)  =  (R,  0)  of  the  boundary 
(My  x  8/(qR2(l  —  u)))  for  second-order  element  simply-supported  plate  solutions  with 
bi-quadratic  mapping  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  1.0  [17]. 
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Figure  18:  Fully-Clamped  circular  plate  model 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.84116 

0.94176 

0.98123 

0.99478 

0.99884 

1.00001 


0.64189 

0.92727 

0.98313 

0.99615 

0.99935 

1.00014 


0.84064 

0.94128 

0.98077 

0.99432 

0.99839 

0.99956 


0.64117 

0.92682 

0.98272 

0.99575 

0.99894 

0.99974 


0.84064 

0.94127 

0.98076 

0.99432 

0.99839 

0.99955 


0.64116 

0.92682 

0.98272 

0.99574 

0.99894 

0.99973 


Table  21:  Center  displacement  (w  x  64 D / (qR^))  for  first-order  element  fully-clamped 
circular  plate  solutions  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  1.0  [17]. 


are  given  by  [17]: 


u,(0,0)=  qR 

(25a) 

M,(  0,0)=  (1  +  0^ 

(25b) 

MX(R,  0)  =  ~ 

(25c) 

Normalizations  are  taken  with  respect  to  these  quantities. 

We  first  study  the  convergence  of  the  center  displacement.  Tables  21  and  22  display 
the  convergence  of  the  normalized  center  displacement  for  the  first-  and  second-order  plate 
elements,  respectively,  obtained  for  the  two  refinement  schemes  with  various  meshes. 

We  next  study  the  convergence  of  the  normalized  bending  moment  about  the  x-axis  at 
point  (x,  y)  =  (R,  0)  of  the  clamped  edge.  We  sample  the  discrete  bending  moment  at  a 


28 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.98752 

0.99984 

1.00043 

1.00046 

1.00046 

1.00046 


1.00406 

1.00066 

1.00042 

1.00040 

1.00040 

1.00040 


0.98707 

0.99939 

0.99998 

1.00000 

1.00000 

1.00000 


1.00366 

1.00027 

1.00002 

1.00001 

1.00000 

1.00000 


0.98706 

0.99938 

0.99997 

0.99999 

0.99997 

0.99995 


1.00365 

1.00027 

1.00002 

1.00000 

1.00000 

1.00000 


Table  22:  Center  displacement  (w  x  64 D/(qR4))  for  second-order  element  fully-clamped 
circular  plate  solutions  with  bi-quadratic  mapping  for  the  three  thickness/length  ratios. 
Reference  thin  plate  limit  solution  is  1.0  [17]. 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.26160 

0.59785 

0.79476 

0.89692 

0.94842 

0.97421 


0.22495 

0.56414 

0.78106 

0.89274 

0.94741 

0.97410 


0.26160 

0.59785 

0.79476 

0.89692 

0.94842 

0.97421 


0.22478 

0.56398 

0.78088 

0.89255 

0.94720 

0.97389 


0.26160 

0.59785 

0.79476 

0.89692 

0.94842 

0.97421 


0.22478 

0.56397 

0.78088 

0.89255 

0.94720 

0.97389 


Table  23:  Bending  moment  about  the  x-axis  at  point  (x,  y )  =  (R,  0)  of  the  clamped  edge 
{Mx  x  (—8 )/(qR2))  for  first-order  element  fully-clamped  plate  solutions  for  the  three 
thickness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 


quadrature  point  lying  closest  to  (x,  y )  =  (R,  0).  Tables  23  and  24  display  the  convergence 
of  the  normalized  bending  moment  obtained  using  the  two  refinement  schemes  for  the 
first-  and  second-order  plate  elements,  respectively.  We  note  that  the  bending  moment 
is  converging  linearly  for  both  the  first-  and  second-order  cases,  with  the  second-order 
case  being  the  more  accurate  in  absolute  value. 

We  also  study  the  convergence  of  the  central  bending  moment  about  the  x-  axis. 
The  obtained  moments  are  presented  in  Tables  25  and  26,  for  the  first-  and  second-order 
cases,  respectively.  We  observe  that,  as  in  the  simply-supported  circular  plate  case  with 
the  polar  discretization  scheme,  the  first-order  approximate  solutions  do  not  converge  to 
the  true  solutions. 

Finally,  we  study  the  convergence  of  the  bending  moment  about  the  direction  form¬ 
ing  45°  with  the  x-axis  at  point  (x,y)  =  (R cos(7t/4),  — i?sin(7r/4))  of  the  clamped  edge 
obtained.  We  only  use  the  SCM  discretization  scheme  in  this  case.  We  note  that  point 
(x,y)  =  (.Rcos(7r/4),  —  i2sin(7r/4))  is  a  singular  point  of  the  SCM  meshes.  The  obtained 
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t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.65090 

0.82572 

0.91282 

0.95641 

0.97820 

0.98910 


0.64522 

0.82675 

0.91369 

0.95681 

0.97846 

0.98931 


0.65090 

0.82572 

0.91282 

0.95641 

0.97820 

0.98910 


0.64509 

0.82660 

0.91352 

0.95663 

0.97827 

0.98912 


0.65090 

0.82572 

0.91282 

0.95641 

0.97819 

0.98908 


0.64509 

0.82659 

0.91351 

0.95663 

0.97826 

0.98912 


Table  24:  Bending  moment  about  the  x-axis  at  point  (x,  y)  =  (R,  0)  of  the  clamped 
edge  (Mx  x  (—8 )/(qR2))  for  second-order  element  fully-clamped  plate  solutions  with 
bi-quadratic  mapping  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit 
solution  is  1.0  [17]. 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.45741 

0.37522 

0.33227 

0.31720 

0.31237 

0.31090 


0.55125 

0.89514 

0.97446 

0.99366 

0.99842 

0.99960 


0.45741 

0.37522 

0.33227 

0.31720 

0.31237 

0.31090 


0.55101 

0.89510 

0.97446 

0.99366 

0.99842 

0.99960 


0.45741 

0.37522 

0.33227 

0.31720 

0.31237 

0.31090 


0.55101 

0.89510 

0.97446 

0.99366 

0.99842 

0.99960 


Table  25:  “Center”  bending  moment  about  the  x-axis  (Mx  x  (16)/(gi?2(l  +  z/)))  for 
first-order  element  fully-clamped  plate  solutions  for  the  three  thickness/length  ratios. 
Reference  thin  plate  limit  solution  is  1.0  [17]. 


t/R  =  0.01  t/R  =  0.001  t/R  =  0.0001 

Mesh  PM  SCM  PM  SCM  PM  SCM 


2x2 
4x4 
8x8 
16  x  16 
32  x  32 
64  x  64 


0.75681 

0.93747 

0.98426 

0.99606 

0.99901 

0.99975 


0.95270 

0.98871 

0.99720 

0.99930 

0.99982 

0.99995 


0.75682 

0.93747 

0.98426 

0.99606 

0.99901 

0.99975 


0.95270 

0.98871 

0.99720 

0.99930 

0.99983 

0.99996 


0.75682 

0.93747 

0.98426 

0.99605 

0.99892 

0.99960 


0.95270 

0.98871 

0.99720 

0.99930 

0.99983 

0.99996 


Table  26:  “Center”  bending  moment  about  the  x-axis  (Mx  x  (16)/(gi?2(l  +  z/)))  for 
second-order  element  fully-clamped  plate  solutions  with  bi-quadratic  mapping  for  the 
three  thickness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 
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Mesh 

o 

o 

o 

o 

o 

I—1 

t/L  =  0.0001 

2x2 

0.65992 

0.66015 

0.66015 

4x4 

0.91912 

0.91945 

0.91945 

8x8 

0.97554 

0.97590 

0.97591 

16  x  16 

0.99264 

0.99308 

0.99309 

32  x  32 

0.99798 

0.99814 

0.99814 

64  x  64 

1.00015 

0.99954 

0.99952 

Table  27:  Bending  moment  about  the  direction  oriented  45°  with  respect  to  the  x-axis 
at  point  (x,y)  =  (Rcos(7t/4),  —  i2sin(7r/4))  of  the  clamped  edge  (Mx  x  (— 8)/(qR2))  for 
first-order  element  fully-clamped  plate  solutions  with  SCM  meshes  for  the  three  thick¬ 
ness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0  [17]. 


Mesh 

t/L  =  0.01 

t/L  =  0.001 

t/L  =  0.0001 

2x2 

0.95362 

0.95391 

0.95391 

4x4 

0.97798 

0.97837 

0.97838 

8x8 

0.99430 

0.99458 

0.99458 

16  x  16 

1.00004 

0.99875 

0.99873 

32  x  32 

1.00163 

0.99982 

0.99970 

64  x  64 

1.00180 

1.00023 

0.99993 

Table  28:  Bending  moment  about  the  direction  oriented  45  °  with  respect  to  the  x-axis 
at  point  (x,y)  =  (i2cos(7r/4),  —  i2sin(7r/4))  of  the  clamped  edge  (Mx  x  (— 8)/(qR2))  for 
second-order  element  fully-clamped  plate  solutions  with  SCM  meshes  and  bi-quadratic 
mapping  for  the  three  thickness/length  ratios.  Reference  thin  plate  limit  solution  is  1.0 
[17]. 


moments  are  displayed  in  Tables  27  and  28,  for  the  first-  and  second-order  cases,  re¬ 
spectively.  As  it  can  be  seen,  the  solutions  are  very  accurate  for  both  the  first-  and 
second-order  cases. 


5  Conclusions 

We  have  generalized  the  twist-Kirchhoff  formulation  of  rectangular  plate  elements  to 
arbitrary  quadrilateral  elements.  The  key  aspect  of  this  generalization  is  the  use  of  Piola 
transformed  Raviart-Thomas  rotation  fields.  This  ensures  that  the  normal  components 
of  rotation  are  continuous  across  element  edges. 

We  have  implemented  the  first  two  members  of  the  twist-Kirchhoff  family  of  elements 
and  studied  convergence  and  accuracy  for  non- rectangular  element  discretizations.  Prob¬ 
lems  studied  included  linearly  elastic  square  and  circular  simply-supported  and  clamped 
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plates,  and  a  simply-supported  rhombic  plate  problem  solved  analytically  by  Morley  [15]. 

In  the  square  and  rhombic  plate  problems  we  employed  asymptotically  parallelogram- 
shape-regular  meshes,  and  satisfactory  convergence  of  all  quantities  considered  was  real¬ 
ized. 

In  the  case  of  circular  plates,  we  studied  polar  meshes,  in  which  the  quadrilateral 
elements  are  degenerated  to  triangles  at  the  center  of  the  plate  creating  a  singularity  in 
the  geometric  mapping,  and  meshes  in  which  a  square  is  mapped  to  a  circle,  for  which 
there  are  four  singular  points  in  the  geometrical  mapping  on  the  boundary  of  the  plate. 
We  were  particularly  curious  about  the  behavior  of  the  bending  moments  at  the  singular 
points.  For  the  second-order  element,  all  quantities,  including  the  bending  moments  at 
the  singularities,  converged  nicely  for  all  meshes.  For  the  first-order  elements,  all  quanti¬ 
ties  converged  nicely  except  for  the  bending  moments  at  the  singularity  at  the  center  of 
the  plate  where  the  quadrilateral  elements  were  degenerated  into  triangles.  We  argued 
on  theoretical  grounds  why  the  first-order  element  cannot  be  degenerated  into  a  triangle 
and  our  numerical  results  for  the  polar  meshes  support  this  conclusion.  A  convergent 
triangular  element  with  exactly  the  same  vertex  displacement  and  mid-edge  normal  ro¬ 
tation  degrees-of- freedom  is  the  famous  Morley  triangle  [16],  which  can  be  recommended 
as  an  alternative  to  degenerating  the  lowest-order  twist-Kirchhoff  quadrilateral. 

Aside  from  this  one  admonition  to  avoid  degenerating  the  lowest-order  elements  to 
triangles,  we  find  the  overall  accuracy  of  the  elements  satisfactory  for  considering  their  use 
for  solving  practical  plate  problems.  The  main  attribute  of  these  elements  is  their  unique 
combination  of  efficiency  and  robustness,  that  is,  they  require  only  r  xr  Gauss  quadrature, 
where  r  is  the  order  of  the  elements,  to  produce  stable  elements  with  no  hourglass  or 
other  singular  modes.  We  believe  this  attribute  is  an  important  and  potentially  deciding 
one  in  explicit  dynamic  analysis  for  which  solution  times  and  storage  requirements  scale 
with  the  number  of  quadrature  points. 
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